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Abstract 

We present new solutions to the strong explosion problem in a 
non-power law density profile. The unperturbed self-similar solu- 
tions discovered by Waxman & Shvarts describe strong Newtonian 
shocks propagating into a cold gas with a density profile falling off 
as r~ w , where uj > 3 (Type-II solutions). The perturbations we con- 
sider are spherically symmetric and log-periodic with respect to the 
radius. While the unperturbed solutions are continuously self-similar, 
the log-periodicity of the density perturbations leads to a discrete self- 
similarity of the perturbations, i.e. the solution repeats itself up to a 
scaling at discrete time intervals. We discuss these solutions and verify 
them against numerical integrations of the time dependent hydrody- 
namic equations. Finally we show that this method can be generalized 
to treat any small, spherically symmetric density perturbation by em- 
ploying Fourier decomposition. 

1 Introduction 

Expanding shock waves are naturally produced by diverse astrophysical phe- 
nomena, such as supernovae, GRBs, stellar winds and more. So far, analyti- 
cal self-similar solutions have been found for several simple cases, of which we 
take special interest in the case of strong spherical shocks propagating into 
a density profile that decays as a power of the radius, p = Kr~^ '. The first 
solutions of this kind to be found, now commonly known as the Sedov- Taylor 
(ST) solution, were given by Sedov, Taylor and Von-Neumann p] [2] [3] for 
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the case u < 3 to describe decelerating shocks. In this paper however we 
shall consider a second class of solutions which was discovered by Waxman 
and Shvarts (WS) for the u > 3 case jl]. In these solutions the shock front 
accelerates because of the rapid decay of the density ahead of the shock, 
causing a part of the flow to be causally disconnected from the inner region 
containing the source of the explosion. Mathematically, the boundary of this 
region appears as a singular point of the hydrodynamic equations somewhere 
between the explosion and the shock, called the sonic point. 

The solutions discussed above, while useful, fall short when describing 
shocks propagating into density profiles that deviate from a simple power 
law decay. This might occur in a variety of astrophysical scenarios, e.g. a 
supernova shock propagating into a modulated stellar wind. For this reason 
it is desirable to generalize as much as possible the external density profile for 
which we can obtain analytic solutions, and this is what we attempt here. It 
should be clarified that while we deal with perturbations we do not perform 
an analysis of stability, but only find solutions corresponding to perturbed 
external conditions. The stability of First type solutions has been studied by 
Ryu & Vishniac [5] and Kushinr et al. [6], and that of second type solution by 
Sari et al. [7], and much of the formalism used for the perturbative analysis 
in this paper is similar these works. 

The plan of the paper is as follows: In section[2]we review the unperturbed 
solutions and the boundary conditions at the shock front and at the sonic 
point. In section [3] we develop the perturbation equations and boundary 
conditions. We then discuss the solutions to these equations and compare 
them to numerical results obtained from a full hydrodynamic treatment of 
the problem. In section H] we discuss a method of generalizing our results to 
accommodate arbitrary small density perturbations, and finally conclude in 
section 

2 The Unperturbed Solution 

We proceed to give a quick review of the unperturbed solutions under con- 
sideration [I]. The physical scenario is the discharge of a large amount of 
energy from a point source at the center of a spherically symmetric distri- 
bution of cold gas. It may be noted that spherical symmetry was chosen 
for its relevance to most astrophysical scenarios, but planar and cylindrical 
geometries may readily be treatd as well. The gas density follows a power 
law behavior such that p(r) = Kr~^ . 
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2.1 The hydrodynamic equations 

We begin with the Euler equations for an ideal fluid with adiabatic index 7 
in spherical symmetry: 

(d t + ud r )p + pr~ 2 d r (r 2 u) = 
p{d t + ud r )u + <9 r (7~Vc 2 ) = 

(d t + ud r )( 1 - 1 c 2 p 1 ^) = 0. (1) 

These equations feature the density p, velocity u and speed of sound c as 
the dependant variables, while the pressure has been eliminated through the 
equation of state p = 7 _1 pc 2 . We will use p or c interchangeably as the 
third dependant variable by merit of convenience. The only relevant scale 
in the problem at late times is the shock radius R(t), and so the self similar 
solutions are given in terms of the variable £ = r/R(t). We define the self- 
similar functions U,C,G and P such that the solutions take this form: 

u(r,t) = MU{£) 
c(r,t)=R£C(0 
p{r,t) = BR~"G(0 

p(r,t)=BR~"R 2 P(0- (2) 

These definitions are supplemented by a scaling law for the shock radius, 
R oc R s . This determines the time dependence of the radius to be 

A(t-t ) a ,5 < 1 
R(t)={ Ae t ' T ,5 = l (3) 

A(t -t) a ,5 <1, 

where a = 1/(1 — 5). The fixing of 5 will be discussed shortly Substituting 
Eqf2]into Eq{T]we obtain two equations for the self similar functions U and 
C: 

dU _A 1 (U,C) 



dlogi A(U, C) 
dC A 2 (U,C) 
dfo^~ A(U,C) ' 



(5) 
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where the functions A, A x and A 2 are 
A = C 2 - (1 - Uf 

A, = U(l - U) ( 1 - U - ^1)- C° Uu 2 I (Q ' - ^ 



a J \ 7 

A 2 = C { (1 - ( ,) ( 1 - r - — ) - l^u r 2(1 - U) + " 1 



a / z \ a 

" C + T^T7 2^ 1 ' (6) 

and an implicit condition for G: 

C- 2 (l - u) x G~*- 1+x £ 3X - 2 = Constant (7) 

where 

A=*±£^. (8) 

Evidently for the solution to pass smoothly through a singular point there 
must exist some £ s where 

A = Ai = A 2 = 0. (9) 
2.2 Boundary conditions 

The boundary conditions at the shock front are determined by the Rankine- 
Hugoniot conditions |8J applied to a strong shock. In terms of the self-similar 
functions these turn out to be 



U{£ = 1) 



2 



7 + 1 



^27(7 ~ 1) 
7 + 1 



G(Z = 1) = (10) 
7 - 1 

The value of 5 has yet to be determined, and to find it we need an additional 
condition. This condition is supplied by the requirement that the solution 
pass smoothly through the sonic point, namely that Eqj9] has a solution. 
Solving this system yields the values of £ s and 5, completing the solution in 
the unperturbed case. In general these values can only be found numerically. 
Type-II solutions with a sonic point have been found and are believed to 
exist for u > ujg(j) > 3. There is a small range of 3 < uj < uj g (j) between 
the first and second type solutions where a third kind of self similar solution 
exists [9J, which are out of the scope of this paper. 
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3 Discretely self similar perturbations 



3.1 The perturbation equations 

We now come to the case of a perturbed density profile. For the perturba- 
tion equations to be tractable we aim at a self similar solution by carefully 
choosing a perturbation whose characteristic wavelength scales like radius. 
Namely, we take the perturbed density profile to be 



p(r) + Sp(r) = Kr~ 



l + o- 



r 



where r$ has dimensions of length and bears only on the phase of the pertur- 
bation. p(r) is the unperturbed density, (3 is the frequency of the perturbation 
and a is a small, real, and dimensionless amplitude. Here and elsewhere we 
take the real part of any complex quantity to be the physically significant 
element. 

The perturbed solution is defined as 

u(r, t) + 5u(r, t) = R£[U(Z) + f(t)5U{£)) 

p(r,t) + Sp(r,t) = KR-"[G(0 + f(t)6G(£)] 

p(r,t) + 5p(r,t) = KR-R 2 [P(0 + f(t)6P(£)}, (12) 

and 

R(t)+5R(t) = R(t)[l + f(t)]. (13) 
To enable separation of variables, the function f(t) is taken to obey 

cl s a fR\ q fR 
d \r J fR 

where q describes the frequency and d the amplitude and phase of the pertur- 
bations behind the shock. Finding the values of these parameters is discussed 
in section 13.21 

Plugging Eq{T2] into EqJT] and taking the first order in f(t) yields the 
self-similar linear equations for the perturbations. This set of equations can 
be written as 

MY' = LY (15) 
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where 





~5G~ 






G£ 




5U 


,M = 





G(u-i)e i 




5P 




.-7^ 






-q-u + 3U + £U' 3G + iG' 

iG{q + 5-1 + 2U + £U') 
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G 








P s 



(t/-l)p' 

p2 J 



(16) 



3.2 Boundary conditions for the perturbations 

At the shock front, the perturbed solution must also obey the Hugoniot 
conditions. From this requirement we find that q = if3 and that the boundary 
conditions for the perturbations are 

5G(£ = 1) = ^±A(d -a)- G'il) 
7-1 

6U(£ = 1) = -^—q - U'{1) 
7 + 1 

5P(£ = 1) = -4r[2(g + 1) - « + <*] - P'(l). (17) 
7 + 1 

In analogy to the unperturbed solution, where the parameter 5 was fixed 
by requiring the solution to pass smoothly through the sonic point, we here 
fix d by the same requirement for the perturbation functions. For a general 
value of d, if we start with EqJTT] and solve back towards decreasing radii 
the solution will diverge at the sonic point. There is, however, one value of 
d where this will not happen, and that is the physical value that we seek. 

We can now see that since the real part of f(t) is periodic, the solution is 
discretely self-similar, i.e. it repeats itself up to a scaling factor in intervals 
of ^ = e~f — 1. While the unperturbed solution and the perturbations in 
their complex form are both self-similar, the physical solution which is the 
real part of their sum is not. 



3.3 The Discretely self-similar solution 

While self similarity simplifies the problem by reducing the partial differential 
equations for the perturbations into ordinary differential equations, still these 
are generally not analytically solvable. Therefore, for each specific set of 
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Figure 1: The self-similar functions 5G, <5£7 and <5P as calculated in a numer- 
ical simulation with 7 = 5/3, uj = 17/4 and /3 = 20 and 40 for the blue and 
green lines respectively. The solid and dashed black lines are the real parts 
of the corresponding solutions of the perturbation equations for (3 = 20 and 
(3 = 40 respectively. 

values for 7, u and (3 we must find numerically the functions SG,SU and 6P 
and the parameter d. This can be straightforwardly done, once d is found, 
by integrating back from the shock towards the sonic point. 

In figure [Tj we present some solutions to the perturbation equations dis- 
cussed above, against a numerical solution of the partial differential hydrody- 
namic equations given by Eq.flTJ, corresponding to the same external density 
profile and adiabatic index. The code we employ uses a second order Go- 
dunov scheme to numerically evolve the hydrodynamic equations. The two 
solutions are identical up to small numerical errors, verifying the validity of 
our method. In figure [5] the solution for 5G is plotted at 20 different times, 
separated by a quarter of the period of the density perturbation. Thus five 
periods of the perturbations are represented, and four different phases within 
each. Clearly the shape of the physical solution as a function of £ changes 
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with time because of the factor R l/3 in the function f(t), and yet it repeats 
itself at discrete periods, making it discretely, rather than continuously, self- 
similar. 

It stands out in figure Q] that the characteristic wavelength of the den- 
sity perturbations is notably shorter than that of the velocity and pressure 
perturbations. This happens because the density supports both traveling 
sound waves and stationary (in the fluid rest frame) fluctuations for which 
the pressure and velocity are not perturbed, while the pressure and velocity 
perturbations must propagate as travelling waves. The dominant compo- 
nent of these perturbations near the shock are left traveling sound waves, 
due to the relative weakness of reflected waves from the hydrodynamic pro- 
files behind the shock. From this argument it follows that the characteristic 

wavelengths are given by — £U + \J 7^) for the pressure and velocity, 
together with ^(1 — £JJ) for density perturbations. 

Finally in figures [3] and H] we look at the parameter d(/3), relating the 
fractional perturbation in the shock position to the fractional perturbation 
in the external density, for several values of u; and 7. It can be seen that while 
the real part of d is roughly constant and of order unity, the imaginary part 
grows approximately linearly with (3. The implication for the perturbations 
is that 5R becomes small (on the order of u/\d\) when f3 is large. This is 
physically sensible since when (3 — > 00 the perturbations oscillate so quickly 
that the shock position, which is the integral of the shock velocity, does 
not respond quickly enough to be significantly affected by the perturbations. 
The other perturbations, 5G, 5U and 5P, are themselves of order d, as can 
be seen from Eq.f fTTl) . and so the actual perturbations to the hydrodynamic 
quantities remain of order a. We can work out the value of the imaginary 
part of d when f3 tends to infinity by considering the short wavelength limit. 
The wavelength of the wave excited by the external density perturbations 
is then very short compared to the scale of variations in the background 
hydrodynamic quantities, and these perturbations can then be approximately 
treated as left traveling waves in a uniform medium. Such waves satisfy 
5u = —5p/(pc) [H], and using equation. (flUj) . (JTTJ) we arrive at the relation 



which for large q gives an excellent approximation to the numerically cal- 
culated value for Im(d), explaining both the magnitude and phase of high 
frequency perturbations. We defer discussion of the small wavenumber limit 
to Appendix [A] 




(18) 
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Figure 2: Discrete Self-similarity: 5G is plotted in four different phases of 
its periodical repetition. The blue, green, red and cyan lines correspond to a 
phase shift of 0,1/4,1/2 and 3/4 period relative to figure [H In each color five 
different periods are drawn, but because of the discrete self-similarity they 
overlap and are almost indistinguishable. 



9 



-3.4 




10 20 30 40 50 10 20 30 40 50 



Figure 3: The real and imaginary parts of d(/3) for several values of uj at 



7 = 1- 




Figure 4: The real and imaginary parts of d(/3) for several values of 7 at 
uj = 5. The slopes of Im(d) computed numerically from this graph are 
(4.830,4.450,4.237) for 7 = (4/3,3/2,5/3), and agree to 3 decimal places 
with those predicted by equation. (|18|) for the respective values of 7. 
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Figure 5: Left: Theoretical and numerical results for a square wave density 
perturbation with (3 = 20 and a = 0.1. Right: The power law external 
density profile perturbed by a square wave, plotted with a = 0.8 to make the 
perturbations more visible. 

4 Arbitrary perturbations 

We have so far discussed the solutions of the perturbation equations for small 
log-periodic perturbations in the external density. However, the linearity of 
the perturbations allows us to construct more general solutions by a Fourier 
decomposition of any periodic spherically symmetric density perturbation. 
We thus treat our basic solutions as a basis for a more general solution 
space. It should be noted that while the basic perturbations are self similar 
(discretely so considering only the real part, but continuously if treated as 
complex functions), a solution that is the sum of solutions with different /3's 
will not be self similar and will have a time dependent profile. This can 
be plainly understood by analogy to the Schroedinger equation in quantum 
mechanics. There, each energy eigenfunction is time independent up to a 
phase, but once different eigenfunctions are combined, their sum becomes 
time dependent owing to the phases of different eigenfunctions changing at 
a different rate. 

We confirm the validity of this method by solving the full nonlinear PDEs 
numerically for a square wave density perturbation. In figure |5] we compare 
this solution to the theoretical solution obtained by summing a large, but 
finite, amount of the Fourier components that constitute the desired exter- 
nal density profile, where each component is calculated using the methods 
described above. This cut off series inevitably creates unphysical oscillations 
(Gibbs phenomenon) in the theoretical solution which may be disregarded. 

This method can in principle be applied even to non-periodic perturba- 
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tions, such as an isolated pulse or step function (density jump). These would 
require a continuous Fourier transform to accomplish which is technically dif- 
ficult when, as is the case here, the perturbations are obtained by numerically 
solving an ODE. 

5 Discussion 

We have laid out a method for solving the strong explosion problem in density 
profiles that deviate from a pure power law radial dependence. The key 
lies in choosing radially log-periodic perturbations which do not introduce 
a new scale into the problem. This leads to self-similar perturbation in the 
hydrodynamic quantities behind the shock, which can be found by solving a 
set of ordinary differential equations. The perturbations are fully self-similar 
when the density perturbation is formally taken to be a complex function, 



Sp/p = o [J^j > but taking the physical real part of the solution makes 

the perturbations, as well as the full solution (the sum of the unperturbed 
solution and the perturbation), only discretely self similar because of the 
periodic nature of the perturbation. We find that the coefficient d connecting 
the amplitude of the perturbations in the shock position with the amplitude 
of the density perturbations has a 0(1) real part and an 0(/3) imaginary 
part, so at short wavelengths, (3 3> 1, this perturbation becomes small and 
is at a relative phase of a quarter wave behind the density perturbation. 

The linearized perturbation treatment naturally ensures that the pertur- 
bations will not depend on a other than by a linear scaling of the amplitude. 
This simplifies the solution of the problem but limits the validity of the 
method to small perturbations. It can be seen in figure El where the solution 
for a specific perturbation is plotted for several different values of a, that 
this limitation becomes pronounced when a « 0.1, although it gives qual- 
itatively correct results even for much higher values. When a is small the 
difference between the simulated and exact values is quadratic with a, which 
is manifested in figure [6] as a linear scaling due to the a normalization of the 
different plots. 

A natural extension of the argument presented here is to cover the rela- 
tivistic regime. The basis for such a study would be self-similar solutions for 
power law density profiles. These were discovered for the ultra relativistic 
limit by Blandford & McKee [10] for first type solutions with u> < 4, and by 
Best & Sari [TT] for second type solutions with u > 5 — a/3/4. An exploration 
of relativistic similarity solutions in various geometries (planar, cylindrical 
and spherical) was given by Sari [TJ]. Additionally Pan & Sari [T3] studied 
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Figure 6: Numerical simulation results for the velocity perturbation at /3 = 
20, 7 = 5/3, tjj = 17/4 and different values of a, divided by a . The difference 
between the lines comes from nonlinear terms that become pronounced as 
cr increases. Small values of a show a good agreement of the numerical 
simulation to the linear approximation of our method. 

the case where a shock traverses a star's interior and then emerges into empty 
space. All of these works are valid starting points for perturbative analyses 
such as the one presented here, and will possibly be pursued in future work. 
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A Appendix: The small wavenumber limit 

In considering the family of solutions presented above, an interesting limiting 
case presents itself in the form of the long wavelength limit, namely when 
the frequency (3 vanishes. In this limit the flow at any stage will converge to 
a self similar form before the density perturbation changes significantly. In 
other words, the solution at any instant should be an unperturbed solution of 
the type discussed in section [21 but with the magnitude K of density profile 
slowly changing with time. The value of K is still not enough to uniquely 
determine the form of the solution, since it is possible that when K changes 
the parameters A and to (see equation. (jBJ)) also change, expressing a change 
of the energy or starting time of the explosion. It can be shown that if we 
take 5 A/ A = a/d and Sto = (while from its definition 5K/K = cr) we 
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obtain a solution of the form 



5G(0 = (d - uj)G(0 - £G'(0 
6U(£) = -££/'(£) 

= (d + 2 -«)P(fl -£/"(£). (19) 

(20) 

This can then be explicitly shown to solve the perturbation equations for 
q = 0. This solution which is zero-order in q still does not allow us to 
determine the value of d at the small (3 limit, as the divergence of the solution 
near the sonic point only appears at first order in q, and to that order the 
equations are not analytically tractable in general. 
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